import java.text.DecimalFormat;

public class Matrix
{
	private double[][] theMatrix;

	private double[][] Q; //For the decompositions LU and QR
	private double[][] R;
	
	private double error;
	
	public Matrix( double [][] matrixIn )
	{
		setTheMatrix(matrixIn);
	}

	public double[][] getTheMatrix()
	{
		return deepClone(theMatrix);
	}

	public void setTheMatrix(double[][] theMatrix)
	{
		this.theMatrix = deepClone(theMatrix);
		Q = null;
		R = null;
	}

	public double getError()
	{
		return error;
	}
	
	public double[][] getQ()
	{
		return Q;
	}

	public double[][] getR()
	{
		return R;
	}
	public double[] getColumn(int n)
	{
	    double [] result = null;
	    if(n < numColumns())
	    {
		result = new double[numRows()];
		for( int i = 0; i < numRows(); i++ )
		{
		    result[i] = this.theMatrix[n][i];
		}
	    }
	    
	    return result;
	}	
	public double[] getRow(int n)
	{
	    double [] result = null;
	    if(n < numRows())
	    {
		result = new double[numColumns()];
		for( int i = 0; i < numColumns(); i++ )
		{
		    result[i] = theMatrix[i][n];
		}
	    }
	    
	    return result;
	}

  public int numRows() {
    return theMatrix.length;
  }

  public int numColumns() {
    return theMatrix.length == 0 ? 0 : theMatrix[0].length;
  }
	
	//--------------- BASIC MATRIX OPERATIONS
  public Matrix matrixMultiply(Matrix m) {
    if (numColumns() != m.numRows())
      throw new IllegalArgumentException("matrix dimension mismatch");
    double [][] result = new double[numRows()][m.numColumns()];
    for (int i = 0; i < result.length; i++)
      for (int j = 0; j < result[i].length; j++) {
        result[i][j] = 0;
        for (int k = 0; k < m.numRows(); k++)
          result[i][j] += theMatrix[i][k] * m.theMatrix[k][j];
      }
    return new Matrix(result);
  }

  public Matrix matrixAdd(Matrix m) {
    if ((numRows() != m.numRows()) || (numColumns() != m.numColumns()))
      throw new IllegalArgumentException("matrix dimension mismatch");
    double [][] result = new double[numRows()][numColumns()];
    for (int i = 0; i < result.length; i++)
      for (int j = 0; j < result[i].length; j++)
        result[i][j] = theMatrix[i][j] + m.theMatrix[i][j];
    return new Matrix(result);
  }

  public Matrix transpose() {
    double [][] matrix = new double[numColumns()][];
    for (int i = 0; i < matrix.length; i++) {
      matrix[i] = new double[numRows()];
      for (int j = 0; j < matrix[i].length; j++)
        matrix[i][j] = theMatrix[j][i];
    }
    return new Matrix(matrix);
  }

   double[] vectorMultiply(double[] vector) {
    double [] result = new double[theMatrix.length];
    for (int i = 0; i < result.length; i++)
      result[i] = dotProduct(theMatrix[i], vector);
    return result;
  }

  public Matrix scalarMultiply(double d) {
    double[][] matrix = new double[theMatrix.length][];
    for (int i = 0; i < matrix.length; i++) {
      matrix[i] = new double[theMatrix[i].length];
      for (int j = 0; j < matrix[i].length; j++)
        matrix[i][j] = theMatrix[i][j] * d;
    }
    return new Matrix(matrix);
  }

  public double maxNorm() {
    double max = 0;
    for (double[] vector : theMatrix)
      for (double d : vector)
        max = Math.max(max, Math.abs(d));
    return max;
  }

  public double determinant() {
    for (double [] row : theMatrix)
      if (numRows() != row.length)
        throw new IllegalArgumentException("matrix is not square");
    double [][] u = matrixLUDecomposition().snd.theMatrix;
    double det = 1;
    for (int i = 0; i < u.length; i++)
      if (i < u[i].length)
        det *= u[i][i];
    return det;
  }

  public Matrix subExclude(int n, int m) {
    double[][] matrix = new double[theMatrix.length-1][];
    for (int i = 0; i < matrix.length; i++) {
      matrix[i] = new double[theMatrix[i].length-1];
      double [] source = theMatrix[i < n ? i : i+1];
      System.arraycopy(source, 0, matrix[i], 0, m);
      System.arraycopy(source, m+1, matrix[i], m, matrix[i].length-m);
    }
    return new Matrix(matrix);
  }

  public double trace() {
    double trace = 0;
    for (int i = 0; i < theMatrix.length; i++)
      if (i < theMatrix[i].length)
        trace += theMatrix[i][i];
    return trace;
  }
  	public static Matrix identityMatrix( int n )
  	{
  	    Matrix result = null;
  	    if( n > 0 )
  	    {
  		double ident[][] = new double[n][n];
  		for( int i = 0; i < n; i++ )
  		{  		
  		    for( int j = 0; j < n; j++ )
  		    {
  			ident[i][j] = (i == j) ? 1.0 : 0.0;
  		    }    
  		}
  		result = new Matrix(ident);
  	    }
  	    
  	    return result;
  	}
	//--------------- PART 1

  public static Matrix hilbertMatrix(int n) {
    double[][] h = new double[n][n];
    for (int i = 0; i < n; i++)
      for (int j = 0; j < n; j++)
        h[i][j] = 1d/(i+j+1);
    return new Matrix(h);
  }

  public Tuple<Matrix, Matrix> matrixLUDecomposition()
  {
    // this is for square matrix only (at least for now)
    if (numRows() != numColumns())
      throw new IllegalArgumentException("not yet implemented");
    Matrix u = new Matrix(theMatrix);
    Matrix l = identityMatrix(numRows());
    for (int i = 0; i < u.theMatrix.length; i++)
      for (int j = i+1; j < u.theMatrix.length; j++) {
        double factor = u.theMatrix[j][i] / u.theMatrix[i][i];
        l.theMatrix[j][i] = factor;
        for (int k = i; k < u.theMatrix[j].length; k++)
          u.theMatrix[j][k] -= factor * u.theMatrix[i][k];
      }
    return new Tuple<Matrix, Matrix>(l, u);
  }
  
  public void householder()
  {
      // this is for square matrix only (at least for now)
      if (this.numRows() != this.numColumns())
        throw new IllegalArgumentException("not yet implemented");
      int num = this.numColumns();
      double x[] = null;
      double y[] = null;
      double h[][] = null;
      Matrix ha = new Matrix(this.getTheMatrix());
      Matrix haList[] = new Matrix[num - 1];
      Matrix tempH;
      for( int i = 0; i < num-1; i++)
      {
	  h = new double[num][num];
          y = new double[num - i];
          y[0] = 1.0;
          x = ha.getColumn(i);
          tempH = householder(chopVector(x,i),y);
          
          for( int j = 0; j < num - tempH.numColumns(); j++)
          {
              for( int k = 0; k < num - tempH.numColumns(); k++)
              {
        	  h[j][k] = (j==k) ? 1 : 0;
              }
          }
          for( int j = num - tempH.numColumns(), l = 0; j < num; j++, l++)
          {
              for( int k = num - tempH.numColumns(), m = 0; k < num; k++, m++)
              {
        	  h[j][k] = tempH.getTheMatrix()[l][m];
              }
          }
          System.out.println("h");
          System.out.println(new Matrix(h).toString());
          haList[i] = new Matrix(h);
          ha = ha.matrixMultiply(new Matrix(h));
          //System.out.println("ha");
          //System.out.println(ha.toString());
      }
      this.R = ha.getTheMatrix();
      Matrix qCalc = haList[haList.length-1];
      for( int i = haList.length-2; i >=0 ; i-- )
      {
	  qCalc = qCalc.matrixMultiply(haList[i]);
      }
      this.Q = qCalc.getTheMatrix(); 
      //return ha;
  }
  
  public static Matrix householder(double [] x, double [] y) {
    x = normalize(x);
    y = normalize(y);
    double [] w = vectorDifference(x,y);
    //System.out.println( "W: " + doubleToString(w));
    Matrix z = outerProduct(w, w);
    z = z.scalarMultiply(-2/lengthSquared(w));
    return identityMatrix(z.numRows()).matrixAdd(z);
  }
  private int lowestNonZeroInVector( double [][] colIn )
  {
	  int result = 0;
	  for( int i = 0; i < colIn[0].length; i++ )
	  {
		  if( Math.abs(colIn[0][i]) > 0.00001  )
		  {
			  result = i;
		  }
	  }
	  
	  return result;
  }
  public void givens()
  {
      // this is for square matrix only (at least for now)
      if (this.numRows() != this.numColumns())
        throw new IllegalArgumentException("not yet implemented");
      Matrix gList[] = new Matrix[(this.numRows()-1)*(this.numColumns()-1)];
      Matrix A;
      int k = 0;
      A = new Matrix(this.getTheMatrix());
      for( int i = 0; i < this.numColumns()-1; i++ )
      {
          for( int j = i+1; j < this.numRows(); j++ )
          {
        	  if( Math.abs(this.theMatrix[j][i]) > 0.00001 )
        	  {
        		  double [][]aMat = A.getTheMatrix();
	        	  double r = Math.sqrt(Math.pow(aMat[i][i],2)+Math.pow(aMat[i][j],2));
	        	  double c = aMat[i][i] / r;
	        	  double s = -aMat[i][j] / r;
	        	  
	        	  System.out.println( r );
	        	  System.out.println( c );
	        	  System.out.println( s + "\n\n" );
	        	  
	        	  double [][]opMat = Matrix.identityMatrix(this.numColumns()).getTheMatrix();
	        	  opMat[i][i] = c;
	        	  opMat[j][i] = -s;
	       	      opMat[i][j] = s;
	       	  	  opMat[j][j] = c;
	       	  	  gList[k] = new Matrix(opMat);
	        	  System.out.println( gList[k].toString() );
	       	  	  A = A.matrixMultiply(gList[k]);
	       	  	  System.out.println( A.toString() );
	       	  	  k++;
        	  }
          }
      }
      k--;
      this.R = A.getTheMatrix();
      A = gList[k].transpose();
      for( int i = k-1; i >=0 ; i-- )
      {
    	 A = A.matrixMultiply(gList[i].transpose());
      }
      this.Q = A.getTheMatrix();
  }
	//--------------- PART 2
	//TODO 
	//--------------- PART 3
	//TODO
	//--------------- PART 4
	//TODO

  public static Matrix randomTwoByTwo() {
    return randomTwoByTwo(-2, 2);
  }

  public static Matrix randomTwoByTwo(double rangeMin, double rangeMax) {
    double[][] matrix = new double[2][2];
    for (double [] vector : matrix)
      randomizeVector(vector, rangeMin, rangeMax);
    return new Matrix(matrix);
  }

  public Tuple<Double,Integer> power(int maxIterations, double threshold) {
    double [] w = new double[numColumns()];
    double [] u = new double[numColumns()];
    double [] ou = new double[numColumns()];
    randomizeVector(w, 0, 1);
    randomizeVector(u, 0, 1);
    double eigen = dotProduct(w, u) / dotProduct(w, ou);
    double oeigen = 0;
    int iterations = 0;
    while (iterations++ < maxIterations) {
      if (Math.abs(eigen - oeigen) < threshold)
        break;
      ou = u;
      u = vectorMultiply(ou);
      oeigen = eigen;
      eigen = dotProduct(w, u) / dotProduct(w, ou);
    }
    return new Tuple<Double,Integer>(eigen, iterations);
  }

  public Tuple<Double,Integer> inversePower(int maxIterations,
                                            double threshold) {
    Tuple<Matrix,Matrix> lu = matrixLUDecomposition();
    double [] w = new double[numColumns()];
    double [] u = new double[numColumns()];
    double [] ou = new double[numColumns()];
    randomizeVector(w, 0, 1);
    randomizeVector(u, 0, 1);
    double eigen = dotProduct(w, u) / dotProduct(w, ou);
    double oeigen = 0;
    int iterations = 0;
    while (iterations++ < maxIterations) {
      if (Math.abs(eigen - oeigen) < threshold)
        break;
      ou = u;
      u = lu.snd.backSubstitute(lu.fst.forwardSubstitute(ou));
      oeigen = eigen;
      eigen = dotProduct(w, u) / dotProduct(w, ou);
    }
    return new Tuple<Double,Integer>(1/eigen, iterations);
  }

  // utility methods
  public static double[][] deepClone(double[][] in) {
		double [][] clone = new double[in.length][];
    for (int i = 0; i < clone.length; i++)
      clone[i] = deepClone(in[i]);
    return clone;
  }

  public static double[] deepClone(double[] in) {
    double [] clone = new double[in.length];
    System.arraycopy(in, 0, clone, 0, clone.length);
    return clone;
  }

  public static void randomizeVector(double [] vector,
                                     double rangeMin, double rangeMax) {
    double rangeSize = rangeMax - rangeMin;
    for (int i = 0; i < vector.length; i++)
      vector[i] = (Math.random() * rangeSize) - rangeMin;
  }

  public static double dotProduct(double[] a, double[] b) {
    if (a.length != b.length)
      throw new IllegalArgumentException("vector dimension mismatch");
    double dp = 0;
    for (int i = 0; i < a.length; i++)
      dp += a[i] * b[i];
    return dp;
  }

  public static Matrix outerProduct(double [] a, double [] b) {
    double [][] op = new double[a.length][b.length];
    for (int i = 0; i < op.length; i++)
      for (int j = 0; j < op[i].length; j++)
        op[i][j] = a[i] * b[j];
    return new Matrix(op);
  }

  public static double[] vectorDifference(double [] a, double [] b) {
    if (a.length != b.length)
      throw new IllegalArgumentException("vector dimension mismatch");
    double [] difference = new double[a.length];
    for (int i = 0; i < difference.length; i++)
      difference[i] = a[i] - b[i];
    return difference;
  }
  public static double[] vectorSum(double [] a, double [] b) {
      if (a.length != b.length)
        throw new IllegalArgumentException("vector dimension mismatch");
      double [] difference = new double[a.length];
      for (int i = 0; i < difference.length; i++)
        difference[i] = a[i] + b[i];
      return difference;
    }
  public static double[] vectorScale(double[] v, double s) {
    double [] result = new double[v.length];
    for (int i = 0; i < result.length; i++)
      result[i] = v[i] * s;
    return result;
  }
  
  public static double[] chopVector(double[] x, int n)
  {
      double [] result = null;
      
      if( n < x.length && x.length > 1 )
      {
	  result = new double[x.length-n];
	  for( int i = n, j = 0;  i < x.length; i++, j++)
	  {
	      result[j] = x[i];
	  }
      }
      
      return result;
  }
  public static String doubleToString( double [] in )
  {
	      String result = "";
	      DecimalFormat df = new DecimalFormat("#.##");
	      for( int i = 0; i < in.length; i++ )
	      { 
			  result += df.format(in[i]) + "          ";

	      }
	    return result;
  }

  public static double lengthSquared(double[] a) {return dotProduct(a,a); }
  public static double length(double[] a) {return Math.sqrt(lengthSquared(a)); }
  public static double [] normalize(double [] a) {return vectorScale(a, 1/length(a)); }

  public double [] forwardSubstitute(double [] vector) {
    // Matrix is assumed to be in lower-triangular form
    if (numRows() != vector.length)
      throw new IllegalArgumentException("solution has wrong dimension");
    vector = deepClone(vector);
    double [] result = new double[numRows()];
    for (int i = 0; i < result.length; i++) {
      result[i] = vector[i] / theMatrix[i][i];
      for (int j = i; j < result.length; j++)
        vector[j] -= result[i] * theMatrix[j][i];
    }
    return result;
  }

  public double [] backSubstitute(double [] vector) {
    // Matrix is assumed to be in upper-triangular form
    if (numRows() != vector.length)
      throw new IllegalArgumentException("solution has wrong dimension");
    vector = deepClone(vector);
    double [] result = new double[numRows()];
    for (int i = result.length; i > 0; i--) {
      result[i-1] = vector[i-1] / theMatrix[i-1][i-1];
      for (int j = i; j > 0; j--)
        vector[j-1] -= result[i-1] * theMatrix[j-1][i-1];
    }
    return result;
  }

  public Tuple<Matrix,Matrix> exciseDiagonal() {
    double[][] d = new double[numRows()][numColumns()];
    double[][] r = deepClone(theMatrix);
    for (int i = 0; i < numRows() && i < numColumns(); i++) {
      d[i][i] = r[i][i];
      r[i][i] = 0;
    }
    return new Tuple<Matrix,Matrix>(new Matrix(d), new Matrix(r));
  }

  public Matrix entrywiseInvertDiagonal() {
    double [][] matrix = deepClone(theMatrix);
    for (int i = 0; i < matrix.length; i++)
      if (i < matrix[i].length)
        matrix[i][i] = 1/matrix[i][i];
    return new Matrix(matrix);
  }

  public Tuple<double[],Integer> jacobiIteration(double [] b, int maxIterations, double threshold) {
    Tuple<Matrix,Matrix> dr = exciseDiagonal();
    dr.fst = dr.fst.entrywiseInvertDiagonal();
    double [] x = new double[b.length];
    double [] ox = new double[b.length];
    randomizeVector(x, 0, 1);
    int iterations = 0;
    while (iterations++ < maxIterations) {
      if (lengthSquared(vectorDifference(x, ox)) < threshold)
        break;
      ox = x;
      x = dr.fst.vectorMultiply(vectorDifference(b, dr.snd.vectorMultiply(ox)));
    }
    return new Tuple<double[],Integer>(x, iterations);
  }

  public Tuple<Matrix,Matrix> splitUpperLower() {
    double [][] l = new double[numRows()][numColumns()];
    double [][] u = new double[numRows()][numColumns()];
    for (int i = 0; i < theMatrix.length; i++)
      for (int j = 0; j < theMatrix[i].length; j++)
        (i < j ? u : l)[i][j] = theMatrix[i][j];
    return new Tuple<Matrix,Matrix>(new Matrix(l), new Matrix(u));
  }

  public Tuple<double[],Integer> gaussSeidelIteration(double [] b, int maxIterations, double threshold) {
    Tuple<Matrix,Matrix> lu = splitUpperLower();
    double [] x = new double[b.length];
    double [] ox = new double[b.length];
    randomizeVector(x, 0, 1);
    int iterations = 0;
    while (iterations++ < maxIterations) {
      if (lengthSquared(vectorDifference(x, ox)) < threshold)
        break;
      ox = x;
      x = lu.fst.forwardSubstitute(vectorDifference(b, lu.snd.vectorMultiply(ox)));
    }
    return new Tuple<double[],Integer>(x, iterations);
  }

  public String toString() {
      String result = "";
      DecimalFormat df = new DecimalFormat("#.##");
      for( int i = 0; i < this.theMatrix.length; i++ )
      {
	      for( int j = 0; j < this.theMatrix[0].length; j++ )
	      {
		  
		  result += df.format(theMatrix[j][i]) + "          ";
	      }
	      result += "\n";
      }
    return result;
  }
}
